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Abstract 

In this paper we consider elliptic and hyperbolic problems in unbounded 
regions. These problems, when one wants to solve them numerically, have the 
difficulty of prescribing boundary conditions at infinity. Computationally, 
one needs a finite region m which to solve these problems. The corresponding 
conditions at infinity imposed on the finite distance boundaries should 
dictate the boundary condition at infinity and be accurate with respect to the 
interior numerical scheme. Such boundary conditions are commonly referred to 
as absorbing boundary conditions. This paper presents a survey and covers our 
own treatment on these boundary conditions for wave-like equations. 


Research was supported in part by the National Aeronautics and Space 
Administration under NASA Grant NAG-1-527 and in part by NASA Contract No. 
NAS 1-17070 while the author was in residence at the Institute for Computer 
Applications in Science and Engineering, NASA Langley Research Center, 
Hampton, VA 23665. 


i 





ABSORBING BOUNDARY CONDITIONS 
FOR EXTERIOR PROBLEMS 

by 

S. [. Ilariharan* 


1. Intr oduct ion 

Many formulations arising from physical nature yield problems in un- 
bounded regions Mathematical formulations of such problems yield govern- 
ing partial differential equations in or near a given domain in such a fashion 
that: i) the equations may be linear but with non-constant coefficients, or 
ii) the equations may be nonlinear, but at large distances essentially be- 
have linearly and with constant coefficients. This note presents a survey 
of the treatment of such problems, when the desired solutions are governed 
by elliptic or hyperbolic partial differential equations These problems are 
called exterior problems and commonly arise in the fields of aerodynamics, 
meteorology, electromagnetic scattering, and atmospheric acoustical wave 
propagation. The mam difficulties with these problems are the boundary 
conditions that need to be prescribed at large distances from the region of 
interest. Usually only an asymptotic behavior is known Such conditions 
may be sufficient to check the well-posedness of the problem however, if 
one wants to compute the solutions of these problems numerically, infinite 
distances need to be truncated to finite distances. The boundary conditions 
imposed on these finite distance boundaries should dictate the behavior at 
infinity and be accurate with the interior numerical scheme. Furthermore, 
the shorter the distances, the more efficient the solutions in terms of com- 
putation time required. This consideration is the essential need in several 
problems, depending on the problem and the kind of computer that is used. 
Typical of such cases are most three dimensional problems. The intention 
here is to present some available techniques to overcome this difficulty, in- 
cluding application to some model problems. 

To begin, let us illustrate some concepts using simple one dimensional 
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model problems. First let us consider an elliptic problem where the field 
equations are reduced wave equations. 


Consider a slab of thickness L with a varying index of refraction n{x) 
Left of this slab, let the media be homogeneous with index of refraction l 
and right of the slab ( x > L ) let the index of refraction be tiq (> l), which 
is a constant Such problems are common in optics as well as geophysical 
waves. Tins situation is illustrated in Figure 6.1 The governing equations 
are as follows. 

a" + k 2 u =0 x < 0 

u" -f- k 2 n 2 (x)a =0 0 < x < L (l.l) 

u" + k 2 riQU =0 x > L 


INCIDENT W/VE 


ARTIFICIAL BOUNDARY 


EXPtikX) 

u" ♦ k 2 u = 0 


REFLECTED 

WAVE 


TRANSMITTED WA/E 


u%k 2 iWu 


0 




Figure 1 One dimensional model 


Then (l.l) needs to be solved in (0, L '\ , ( L ' > L ) say with boundary condi- 
tions 

0|»(0)=J (12) 

IJjii(L') = 0. (1.3) 

u, u' continuous on interfaces, (1 |) 

where D\ and Bt are boundary operators and g are given are to be chosen 
according to the physics of the problem. For example, if there is a unit 



amplitude wave traveling from — oo incident on the slab, then for x < 0 the 
solution may be written as 

u(x) = e lkx + li{k)e - ,kl (l 5) 

where R(k) is the reflection coefficient and measures the part of the wave 
which is reflected from the boundary x = 0 Eliminating R(k) from (1.5) 
yields 

u'(0) + iku{0) = 2ik (16) 

which is called an inllow condition. Thus in view of (1.2) B[ = d/dr 4- 
ik and y = 2ik. Now we can do a similar calculation to obtain a boundary 
condition at U For x > L all the waves transmit and do not reflect back as 
there are no other boundaries for x > L. In this region the solution may be 
written as 

u(x) = T(k)e ,knnI (1.7) 

where T(k) is the transmission coefficient which measures the transmitted 
part of the wave. Eliminating T(k) in (1.7) we obtain 

a'(L') — iknpu(L') = 0. (1.8) 

Again comparing with (1.3) we see that the operator S 2 is 

d 

Bi = ikn 0 . (19) 

ax 

The boundary condition (18) is the desired absorbing boundary condition 
for this problem, which is exact. 

VVe note that this condition could have been applied exactly at x = L. 
and this, as we will see. does not always hold in higher dimensions. There is 
another lechnital difficulty in this type of problem. One can find that even 
if n 2 (x) = constant there is a countable set of wave numbers {k n } called 
the interior resonant values for which the continuous problem does not have 
solutions. When one deals with higher dimensional problems, it is not easy 
to calculate such frequencies when one has a general shaped region where 
the solution is sought. 
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Figure 2 One dimensional model 


As an example of hyperbolic problems and associated boundary condi- 
tions, let us consider again a one dimensional model In particular let us 
consider the following problem (see figure 2): 


I 

-57-TttM = 0 < x < L 

c 2 (x) 

l 

rjlLll Hxl X ^ Lj 

C 0 


(1 10) 

a, u x , u t continuous on x = L 


(l II) 

u{ 0 ,t) = f{t) 


( 1 12) 

Bu{L',t)= 0 (£' > L). 


( l 13) 

This problem tells us that the wave is propagating from a source at x = 0 
and the wave speed changes in 0 < x < L and for x > L the waves do not 
reflect, in this region. Equation (1. 13) dictates a no reflection condition in 
this region and B is an operator which is to be determined to fit this physical 
nature. For x > L the outgoing wave solution ran be written as 

u(x. t) = o(x — Cot) 


(III) 

which yields 

<’o u x + ut = 0 • 


(1 15) 

Thus the absorbing condition to be imposed at x = i 

' is (1.15) and 


V-w 

III 

+ 

Sri 


(1 16) 
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We see from both one dimensional examples that the absorbing bound- 
ary conditions are exact and easy to obtain. The difficulties arise in higher 
dimensions as we will see. Nevertheless the first problem we described here 
serves as a model for inverse problems that can be investigated numerically 
to calculate n 2 (x), if the no reflection condition B2U = 0 is given. This is 
described in Dunn and llariharan^l 

The approach taken is as follows- In section two we describe a local 
boundary condition procedure for exterior elliptic problems. In section three 
we describe a nonlocal boundary condition for similar problems In section 
four we describe a new method that combines the previous two procedures 
but is more accurate than either. Finally, in section five we describe 
the methods known for tune dependent problems and describe their use for 
a practical problem. 

2 . Local B ound ar y Conditions Procedure 

We pointed out in the last section that the absorbing boundary condi- 
tions are exact in one dimension and do not hold in higher dimensions We 
want to investigate this statement a little further and provide a treatment for 
this, that of Bayliss, Grunzburger,and TurkelJ 21 For illustrative purposes let 
us consider the following problem that commonly arises in two dimensional 
acoustic scattering. 

Let (2 be a simply connected, bounded domain with boundary T and its 
exterior fi f . Then the problem is (see figure 3 ): 

A u 4- k 2 a = 0 in f2 (21) 

du „ . 

— = 17 on T (2 2 

cm 

u satisfies a radiation condition (2 3 ). 


Briefly this problem has a physical meaning that there is a wave incident on 
fl whose boundary V is a perfect reflector and the reflected waves all decay at 
infinity, u measures the scattered part of the wave and g is the contribution 
from the incident wave. Regardless of the boundary condition on F, the issue 
here is the radiation condition. Note that the radiation condition is another 


RADIATED WAVE 



Figure 3 Two dimensional scattering problem 


name for the absorbing condition It is known that at large distances the 
solution of (2. L) behaves like 

M*) , HO) , 

r r 2 

( 2 . 1 ) 

The first part of this solution dictates outgoing waves and the second part 
dictates incoming waves, which can easily be seen from the signs e~ lkr Thus 
in order for the radiation condition to be satisfied we must pick only the first 
part, i.e. , 


,tkr 


u ~ —j=r ^a o (0) + 


ai{0) a 2 (0 ) 

' 1 


+ 


, - ikr 




bom 


a ~ 



^ao(0) + 


a l {0) 



+ • 


(2 5) 


Comparing this to the one dimensional case in (1.7), shows that instead of 
a single transmission coefficient we have many coefficients a,,i = 0. 1 ? 
Thus the one dimensional way of eliminating the transmission coefficient 
does not exactly work in this situation This shows the difficulty in higher 
dimensions. However, if we settle for loss of accuracy, in particular eliminat- 
ing ao(0),we have 


,ikr 


a r - iku = - K - (ciq(0) + - 

2 r » r 


«l (0) , a 2 {0) 


+ - j- + 


Thus we see that 





(2 6 ) 
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r ji (ii r — iku) = 0(1 /r) 


(2.7) 



or simply ra(u r — iku) = G which is Sommerfeld’s radiation condition 

in two dimensions. 


Taking a closer look at the right hand side of (2.6), in particular, the 
first expression, we see that it is nothing more than —a/2 r. Thus we have 


ll r 


-iku+ -u = 0 
2 r 



( 2 . 8 ) 


Comparing (2 7) and (2 8) shows the error introduced in the boundary con- 
dition at large distances is reduced by a factor of l/r This was first studied 
in connection with numerical implementation for time dependent problems 
by Bayliss and Turkell 3 ! and also by Engquist and Majdat 7 ! For time har- 
monic cases with time dependence of the form e ~ lkt , both of their boundary 
conditions reduce to the form (2.8). Bayhss, Gunzburger,and Turket^ gen- 
eralized these conditions. They defined the boundary operator on the left 
side of (2.8) by 

01 = — -ik + l/2r. (2.9) 

dr 

When (2.8) is explicitly written it has the form 


B i u = 



flt(fl) 2 a 2 {9) 

r 2 + r 3 



( 2 . 10 ) 


They observed that now the second coefficient a\($) can also be eliminated 
If one does that by setting v = B\.a it is readily verified that 

-ikv+ % = 0(l/r 9/2 ) 


or 


/J'2 u — f — — i k + 



(2 U) 


This process can be repeated to generalize 


(2 12 ) 


B m u = 0(l/r 2m+l/2 ) 



where 


D r 


= n 


or r , 


(2 13) 


Thus in order to implement this approach the problem can be considered in 
a truncated region (figure ( f)) as follows: 

Seek tv such that 


A iv + k 2 w = 0 in Hr 

(2.11) 

dw 

(2.13) 

^ = go " 1 

B rn w = 0 on Too 

(2 16) 



Figure 4 Computational domain 


Where is a circle where the approximate boundary conditions are 
to be imposed. 

In the region Ut one can use a finite element formulation of the problem 
similar to the one used in the next section to seek the solution But there 
are some difficulties. 


First, the governing equation is a second order differential equation 
Even the second order boundary condition involves second order derivatives 
in the radial direction. And all the higher order conditions (2. l.'J) have higher 
order radial operators. This can be overcome by appealing to the differential 
equation itself. Rewriting the equation in cylindrical coordinates 


( A 4- k 2 ) w 


0 2 w 1 dw l d 2 w 2 


8 


(2.17) 



From (2.17) we see that second order radial operators can be translated into 
single first order radial derivatives plus derivatives in the tangential direction. 

The second difficulty is the following. In comparing (2.1) - (2.3) and 
(2. 14) - (2. 16) we see that there is an a prim error introduced (irrespective 
of the numerical scheme used) of the order 0( l/R im M / 2 ), where R is the 
distance of I'oo from the origin. This of course depends on the order of the 
operator too. The question is what can we say about || u — w || in a suitable 
norm The answer is not available in two dimensions. However, Bayliss, 
C»unzburgor,and Turkel were able to prove for corresponding operators B m 
in three dimensions for rn = l and rri — 2 the following theorem. Let 
w = a — then 


II ® II ( r) = C/r m + l , (2.18) 

where C depends on k and Too and the norm is defined on the surface as 

II w ||( r) = / M ’ dA . 

This theorem tells that the error on the artificial boundary in that given norm 
is inversely proportional to l/r m + l and this error remains the same in the 
interior boundar> too. The key point here is that this bound is dominated by 
C which depends on the wave number k. For three dimensions as k — * 0 this 
estimate is still valid. However, in two dimensions, even though it is possible 
to get such an estimate as the wave number becomes smaller, this constant 
can grow larger. In fact, in two dimensions there is a logarithmic branch 
point as k — * 0 and the constant C can grow very large. For this reason, for 
low frequency cases the problem must be examined very carefully. Such a 
treatment is technical and it is available in Hariharan and MacCamy 121 

These difficulties are pertinent to the reduced wave equations. However, 
other strongly elliptic cases can be handled without difficulties. For example 
if we consider Laplace's equation Au = 0 ? then the solution at infinity either 
has to be bounded or behaves logarithmically with a given behavior. To 
apply the above process let us consider the following problem: 

Au = () in fi + (2.19) 
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u = g on T 

u ~ log r as t — \J x 2 -\- y 2 


oo. 


( 2 . 20 ) 

( 2 . 21 ) 


Solutions lor this problem can be written as 


OO 

u(r,0) — logr+ ^ ~ cosnO . (2.22) 

n-0 r 


This can be used to obtain boundary conditions on the artificial boundary 
Too. Set v = u — logr — « 0 . Eliminating a t we see that. 


£ + --o(±) 

dr r r J 


Similarly eliminating ci 2 we obtain 


(2-23) 



(2 24) 


This process can be repeated to obtain higher order conditions. The constant 
ao can be obtained by averaging the solution on Too. In section 4 we provide 
an alternative treatment to resolve the difficulties. In fact we implemented 
the conditions (2.2*5) and (2 2 1) along with ours to compare the efficiency of 
our procedure. 


3* Nonloc al Uoundarv Co n ditions Procedure 


This procedure is a little more difficult to handle than the previous 
one. Hut, it does not introduce an a priori error due to placing the arti- 
ficial boundary at finite distances. There are two different versions of this 
procedure. The first one originated in the work of Fix and Marin 191 which 
was done for situations of wave guides and was made .general in two dimen- 
sions by MacCamy and Marird 15 ! The second version is due to Johnson and 
Nedelec' 131 , which was done independently but has similarities in the ap- 
proach. Extension to three dimensions was done by Aziz and Kellogg 1 All 
these procedures are done in view of implementing finite element methods 
As a result, a reasonable amoiint of analysis is available for this method. It is 
almost impossible to summarize all of it here. We choose the method of Mac- 
Camy and Marin and describe how the nonlocal conditions are treated. The 

10 




Figure 5 Two dimensional interface problem 


attractive feature of this method is that it is appropriate to interface prob- 
lems analogous to the one dimensional model that we described in section 
one. Let us present the model problem in two dimensions. Again this is a 
scattering type of a problem, but this time the scatterer is an inhomogeneous 
penetrable body. The problem is as follows. 

A n + k 2 n 2 (x)u = f in 11 (3.1) 

A u -I- k 2 u = 0 in fl + (3.2) 


u 

du~ 

dn 



(3.3) 

(3.4) 


u — u, satisfies Somtnerfeid’s radiation at infinity. (3.5) 

Equations (3.1) and (3.2) say that the index of refraction of the media in 11 is 
n 2 (x),x = (x,//) and outside is just a constant. In (3.5) u,(x) is a prescribed 
incident field which satisfies (3.2). (3.3) and (3.1) are continuity conditions 
of the solution on V In many applications they may not be continuous, but 
rather may have jumps. These jumps can be treated without much difficulty 
To describe the procedure let us extend n 2 (x) in Z 2 so that 


^(x.y) 


n 2 (x.y) in 11 
t in fl + 


and extend / in a similar manner, so that 



in 11 

in 11+ * 


1 1 



The governing equations will thus be 


A u 4- k 2 fj.u = g in Z‘ 


(3.6) 


To do numerical calculations let us truncete the infinite region by a circle 
Too into a finite one as depicted in Figure 6. Denote the truncated region 
by Hr- 



Figure 6 Computational domain 


A standard way of handling this problem is to use variational methods 
Suppose a is a solution and u is an arbitrary differentiable function in Dj 
If we multiply (3.6) by u, integrate over Qt and use integration by parts we 
find , 

- Vu • Vu + k 2 I nuv + / v— = / gv. (3 7) 

Jq t Jqt J r* " n Jq t 

This variational form is the basis of a numerical method. But, in order 
to implement this numerically, one should have enough information about 
Ou/dn on Toq. This in fact requires the knowledge of the solution exterior 
to Too in the region /loo, recalling that in this region a satisfies 

Aa + k 2 a = 0 . 


This suggests that one may seek an integral representation of the solu- 
tions. In particular let us consider 


it(x) = / a{ y) C k {x,y)ds y + u,(x ) 

J r „ 


G l 


OO 


(3.8) 
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where Gk is the free space Green’s function and has the form 


Cjt(x,y) = - l - //q 1) (A: i x - y |). (3.9) 

Tiie representation (3 8) is known as the simple layer representation, winch 
satislies the Helmholtz equation and satisfies the radiation condition at in- 
finity. These two facts arc easily verified by applying the Helmholtz operator 
to (3.8) and by realizing the asymptotic behavior of (3.9) for large values of 
| x | which yields the form (2.3) In order to calculate the normal deriva- 
tive of u on I'oo, (3 8) together with the standard jump relation m potential 
theory yields 

du , , l , , /’ . , d _ , , , duAx.) „ 

<M {k)= 2 a(X)+ Jr CT(y) ^ G ‘ (X ’ y) ^ + "an X sr ~- (310) 

Thus Ou/On ran be calculated once a is known. One can use (3.8) to obtain 
a singular integral equation of the first kind to determine a when x is on 
Too- 

a s = u - u, = / cr(y)Gfc(x, y)ds y , x E V ^ (3.11) 

J r» 

Let us denote this equation in an operator form 

u,(x) = G k o <r(x), xeToc • (3.12) 

Then a can be formally inverted to obtain its value of the form 

cr = C fc ' l ou s . (3 13) 

But there is a technical difficulty here. G exists provided — k 2 is not an 
eigenvalue of the operator A with zero Dirichlet condition on Too- This result 
is familiar in the diffraction theory and it is analogous to the statement that 
we made for the corresponding one dimensional problem. That is to say. 
such values of A- will be the interior resonant values of Hr- But there is 
a way to tre.it this problem. UrselO 181 proposed to seek solutions in the 
exterior not only by a simple layer operator as in (3.8), but combined with 
a double layer representation. This slightly complicates the explanation of 
the present method, but nevertheless can be done. For the moment we shall 



assume k is not a resonant value so that (3.13) holds. Then (3.10) takes the 
form (for x G Too): 


^(*)= ' ° u ’( x ) + f O k { oH 3 (y)^(x,y)ds y +~^-.(3.l[) 

** F OO 


The riglit hand side of (3. IT) is another integral operator acting on u, so 
that (3.1 I) can be written (for x G Too) as. 

dn 

< ')u /On = 7V(u) - T k {u t ) - -r 1 . (3.15) 

an 

Thus, the normal derivative of u is a functional of u on A closer look 
at (3.M) shows that to calculate du/dn at each point on Too we need the 
knowledge of a on the entire boundary Too- Thus the boundary condition is 
nonlocal. Properties of the operator Tfc are discussed in detail in reference^ l5 l 
In particular it is a bounded linear operator. 


Returning to the variational formulation, (3.7) together with (3.15) one 
has 

a(n. v) = - Vu • Vu + k 2 I /.iuv+ v 7\(u) = 

J n t* J o t n xi 

I gv+ [ [Tki'h) + ~ )y = F(t') • (3.16) 

Jn T on 

To provide a brief numerical implementation of (3.16) we seek approximate 
solutions u h such that 

ci(ii h ,v h ) = F(v h ) (3.17) 

for all i'h G S h . 

where u h and e/t are in a finite dimensional subspace S h of an infinite di- 
mensional space S whole u is sought. Then (3 17) is made equivalent to a 


matrix problem by selecting a basis {pi,p 2 > 
can be approximated by 

p\ } for S h . This sajs u 

V 


u h = 

(3.18) 

j = i 


which satisfies 


a(u h .p,) = F(p t ), i = l.* 

• N. (3.19) 

11 




Then (3. 19) is the matrix problem 


k q = f 


(3.20) 


where q = (< 71,1721 • • -,<h\f) T is the vector of weights in (3.20) and f is the 

source term. 

k l„ = <1 {‘Pji<Pi) = - • Vyp t + k 2 / ntp J p l + <PiT k {Pj). 

J n t * n t •' To © 

(» 21 ) 

We can now describe how the nonlocal boundary condition (3.15) is used. 
We see fiotn (3 21) that calculating (3.15) is equivalent to computing the 
integrals 

I <p t T k {<pj)ds (3.22) 

J Tco 

for the basis functions 'pi,'P 2 ,- * ”Pn of the approximation space S h . This 
computation is carried out in a straightforward manner using (3.12) and 
(3. 1 1). First solve for a,, i = !,•••• N. 


J (T l {y)G k {x,y)ds y = *?,(x) , x E t\ 


(3 23) 


and compute Tjt(p,)(x) (for x E f^) from 


\ f Q 

Tk{Pt){x) = &i{y) ic{x,y)ds y . (3.2-1) 

There are elfective procedures to implement (3.23) and (3.24) in two 
dimensions which can be found in MacCamy and Marin 1 and also in 
conjunction with an integral equation treatment to this problem when n 2 (x) 
is a constant found in Kariharan and MacCamy 1 121 

Some final remarks are needed. We assumed u and du/dn are continu- 
ous on the boundary f. If they have jump discontinuities of the form 


/ On \ ( du\ + 

“ 1 ( J7,J • 


(3.23) 



The variational form should be modified. These modifications can be found 
in reference! 15 l and in Bielak, \lacCamy,and McGhee^l together with nu- 
merical implementation. 

4. Infinite ord er radiation condition procedure 

In this section we provide a description of a new method due to Canuto, 
Hariharan,and bustmanI 5 L From the discussions of local boundary conditions 
procedure, we saw that errors introduced in computations are twofold The 
first one is the error which varies according to the farfield distance and the 
order of the boundary condition used. The second one is due to implementa- 
tion of finite element method. In the nonlocal boundary condition procedure 
we see that error is essentially only due to the finite element method and 
the solution of the integral equation on the boundary Too. No other error is 
introduced. The local conditions can be improved by either increasing the 
distance of the artificial boundary or increasing the order of the operator 
B rn so that the dominating error will be essentially due to the finite element 
implementation. Error estimates, in such procedures are usually at the best 
of 0(h 2 ) for this type of problem. A question which then can be asked is 
whether we can improve the accuracy of the solutions by a method which 
does not have an a priori error due to placing the artificial boundary 
and at the same time achieve error in computation 0(h M ) , where M is the 
number of elements used. It is possible, if one uses spectral methods, 
and a brief theoretical discussion of obtaining such accuracy is found 
in Lustman [14] for simple one-dimensional problems. The method we will 
describe here does not have any theoretical error estimates yet. At this 
point we only have numerical evidences. Theoretical error estimates are 
available only in one dimenison, and extension to higher dimensions is 
still much in need of work. The procedure we are going to describe here 
works for a general second-order elliptic problem, exterior to a given 
domain in two dimensions. The structure is as follows: 


Lu = f 
u = g 
Boo 1 * = 0 


in 


on f 


as 



oo , 


(I I) 

( I 

( I ■•!) 


whore and l are same a* that described in section 2. B ^ is a boundary 
IG 



operator at oo. For illustrative purposes and to appraise numerical results 
let us consider the problem given through (2.19) - (2.21), which we repeat 
here. 


Au = 0 in n + (4.1) 

u = g on T (4.5) 

u ~ log r as r — \J x 1 + >j 2 — ► oo (4-6) 

Again the* goal here is treat the condition (4.G) numerically Solution in the 

farfield may be sought through separation of variables of the form 

u(r,ip) = log r + ~^\ e ' kp (4- 7 ) 

k 


where (r, p) are the polar coordinates m the plane. Note that the right hand 
side of ( 1.8) satisfies the radiation condition (4.6) at infinity The coefficients 
ak are unknown. The approach here consists of expressing each coefficient 
a*, as a functional of u, rather than eliminating a finite number of them using 
differential operations described in section two. From (4.7) we see that for 
any r > 0,a/t/r! fc l is the A:- th Fourier coefficient of the periodic function of 
p ► u(r, <p). Thus we can invert the coefficient to obtain 



If we differentiate (1.7) with respect to r 


“r{r,p) 


l 

r 


‘-E i*i 

k 


-!l e ik -p 
r !fc I 


( 1 . 9 ) 


and use ( 1.8), we obtain an integro-differential relation on circle of radius r. 
as follows: 


>lr{r,p) 


I- l r T. 1*1 t**-” <‘(r,e)M 
J 0 . 


( 1 10 ) 


«r 


- K o u 
r 


or 


r 


(bit) 
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where K ou is the convolution of u with singular kernel 


l 00 

K(rj) = - ) k cos krj . 

7T T — 4 


(4.1?) 


fc=i 


Again we consider a truncated region fir with an artificial boundary 
Poo where the condition corresponding to that of (4.11) will be imposed. 
Note at this [joint that (111) is similar to the condition (3 15) except it 
is imposed on a circle of radius r = R^. It is not a necessary restriction 
that Poo should be a circle. But the restriction guarantees accuracy, if the 
approximate solution is a trigonometric polynomial on Too, which is the case 
if a spectral Fourier method is used in the angular direction Suppose the 
approximate solution on l’oo to be u iV , which is a trigonometric polynomial 
of degree V 

* 

=£ l‘l (fi~) ' ,tV ■ (J‘0 

|*|<JV 

The asterisk in the summation indicates periodicity in the <p direction (i e., 
iVy(/?oo) = u^/v(fioo))- Then (1.16) says that the integral operator K pro- 
duces a new polynomial of degree N, whose Fourier coefficients are obtained 
from those of u N (R 00 , •) by multiplication by the modulus of the wave num- 
ber | k | . If it* V is known on F^ through its values a N (/?<*,, £> ; ) at equally 
spaced points jV/A', j =(),!.•• •• ,2 N — l , then K o u N can be computed 
at the same nodes exactly and efficiently, first, by Fourier transforming the 
values of a v to get its coefficients, then multiplying by the modulus of the 
wave numbers and finally by Fourier transforming back to get the point 
values of K o u . This takes order of A r log 2 A r operations if one uses fast 
Fourier transforms, which are always used in this type of calculations. Thus 
the spectral solution is required to satisfy the radiation condition 

u r v = - 1 - [l - K o a"] on T*, . (1 15) 

/too 

Unlike the family of conditions described in section 2, this is precisely the 
same boundary condition satisfied by the exact solution except for the trun- 
cation error which comes from using a finite number of modes. No other error 
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PHY5ICAL DOMAIN COMPUTATIONAL DOMAIN 

(r,0) (s,0) 

Figure 7 Mapping of physical domain to computational domain 


is introduced, and this is why we call this form the infinite order radiation 
condition. 

The actual implementation in a simplified form is as follows: 

Suppose the boundary T is polar representable as r = R{6). We use 
Chcbyshev polynomial approximation in the radial direction and Fourier de- 
composition in the angular direction. For this we need to make a coordinate 
transform as depicted in (igure 7. In particular the transformation is 




So that at s = — l , r(— 1, 0) = R{9). Demanding r( 1, 9) = R^, we determine 
the stretching parameter a(9) as 

«(») = .) log (R„/R(0)) . (1.19) 

This transformation in turn changes the Laplace equation into the Form 

fjU = An = au s , + bu s g + cugg + du s = 0 (4-20) 

where a, 6, c and d are functions of s and 0 obtained from chain rule. Seek 
approximate solutions of the form 

M 

£ £ h m kT m (s)e M (1.21) 

m=0 |fc|<AT 

where T m is the tilth Chebyshev polynomia! of the first kind, defined by 

T m { cos 9) = cos (m9) (4.22) 

and u‘ v is determined by at the ( M + 1) X 2 N Chebysliev- Fourier nodes 

(s t ,9j ) = (cos nr/ M, jtt/ N) ( 1 23) 

i = 0, • • • • M 
j = 0, • • • 2iV - 1 . 

At i = M the given Dirichlet condition should be imposed. At i = 0 ,du/ds 
obtained from Da/Or should be updated. Placing the inhomogeneous terms 
arising from the boundary conditions in a forcing term /, let us rewrite the 
final spectral operator as 

(L, p u‘ v - /) (s.,0,) = 0. (4.21) 

The matrix formed bv L, p a v is large and does not have sparse struc- 
tures. Tims in order to solve ( 1.21), an iterative procedure is desirable. We 
outline only the bruT idea of implementing this procedure. Dropping the 
superscript N in (1.21) the method is as follows: 


(1.25) 


20 


Uri t- 1 — iiti "b r* n ( L s p u n y) 



where a n is chosen so that the 9% norm of the residual 


^nf 1 — f L 3p u n (4. 20) 

is minimum. This gives 

( r n> f J ip r n) 

a n 7~r r T 

V ^isp^ni vipTn) 

where ( , ) denotes an 9 j inner product 

As a sample comparison we can" generate a situation where the exact 
solution is 

a(x, y) = log ((x - -7) 2 + y 2 }. (4.28) 

For the geometry of a circle (r = R{9) = l) and for a 33 x 32 grid, £2 
errors are listed for different values of R^. Solutions are compared again: t 
Bayliss, Gunzburger,and Turkel’s procedure (BGT) with their first order 
(FO) and second order (SO) boundary conditions given in section two (Equa- 
tions 2.23). 10 01 IL denotes the infinite order radiation condition of the 
present method. (Canuto, llariharan and Lustman). 


Roc 

FO BGT 

SO BGT 

10 CHL 

1.2 

.10 x 10“ 2 

.031 x 10~ 2 

.00037 x 10 

3 

.15 x 10~ 3 

.025 x 10~ 3 

01 x 10" 5 

5 

.14 x 10~ 4 

.16 x 10~ 4 

11 x 10~ 4 


TABLE 1 



We see from this table that in all cases the infinite order radiation 
condition is superior, especially when the artificial boundary is near the body, 
which is desired anyway. As R ^ increases the first order and second order 
conditions become better and comparable to the infinite order condition. 
This is because the grid size is fixed If the grid size is increased, the infinite 
order condition will improve. Further illustrations can again be found in 
reference ! 15 L 
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Note further that this procedure can be easily extended to other types of 
problems. If we consider the Helmholtz equation A u+k 2 u = 0 with radiation 
condition at inlinity, we can write the oul going solutions as follows: 

OO 

^ (r,0) = £a n f/<'>(< : r)e‘" ^, (4.30) 

n— 0 


where are the Hankel functions of the first kind and order n. 
If we take the radial derivative of (4.30) we have 


= -SI) 


r»— 0 



Again we note that a n fln\kr) is the nth Fourier coefficient of u(r, 

•), which 

can be inverted to give 



“nUi'Hkr) = £ j 

r2n 


/ u(r,<p)e~ lTt{p d<p . 

0 

(4.32) 

Substitution of this in (1.31) gives 



du y, ktl { n v {kr) 

Ul'Hkr) 

f u(rpp) e tn ( e ~ d<p 

Jo 

( 1.33) 

or 

du 

dr 

K o u . 

(4.34) 


This is similar to what we had in equation (4 LI) except that a finite number 
of values of llj^ and its derivatives must be calculated. To do this there are 
well known recurrence relations and numerical evaluations available in the 
literature. 

5. Bounda r y Co ndi tion s f or Ti m e Dependent Proble ms 

This has been the most difficult problem to handle numerically. It is 
still an open question if a nonlocal condition that is suitable for numerical 
calculations can be obtained. We saw in section 1 that it is relatively easy 



to obtain an absorbing condition for a simple wave equation in one dimen- 
sion. We proceed in the same fashion to obtain boundary conditions in two 
dimensions. This process results in the procedure introduced by Engquist 
and Majdal 7 ^ 8 !. Mefore we proceed to discuss the two dimensional case, note 
that in three dimensions D’Alambert’s solution holds. If we consider 

^2^tt = ^XX Uyy -j- !J ZZ (5.1) 

then the general solution for spherically radiating and incoming waves can 
he written as 


u(r,t) = 1 f(r - ct) + - g{r + ct)- (5.2) 

r r 

This representation is well known and may be found in Morse and Feshback’ ’ 6l 
The first part indicates outgoing waves. If one wants to impose a radiation 
condition, this part may be used to represent the solution as follows: 

u(r,t) = - f(r - ct). (5.3) 

T 

If we assume the sound speed is normalized to one, then it is easy to check 
from (5.3) that 

II T + U( + _ =0 (5.4) 

r 

which is local, both in space and in time, but provides an exact description 
of spherically radiating waves. 

Thus in three dimensions a radiation condition similar to what we dis- 
cussed in sections 2 through 5, in the time domain,will become exact if one 
imposes (5.1) at sulliciently far distances. In two as well as in three dimen- 
sions analogous to (3.2) one wants a cylindrically or spherically radiating 
waves at finite distances and that causes problems, as we will see a little 
later. 

Let us begin the discussion of two dimensional wave equations. Consider 
a wave traveling from left incident on the boundary x — L[L > 0) without 
reflecting (see figure 8) and governed by 


lift — ll’XX -f" Uyy • 


(3.5) 
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(5.5) with respect, to t and y and call the corresponding dual variables r and 
This gives 

u xx + (r 2 - <; 2 )u = 0 . (5 9) 

Equation (5.9) is a reduced wave equation in the x direction, whose solutions 
may be written as 


u(x;r,c) = .'l(r,f)e'v' ; ’-' ri + fl(r, f ) (5.10) 

^ .+ 

Arrows in (5.10) indicate the waves moving in the right and left directions 
In order to impose an outgoing solution we see that we must choose 

u(z;r,f)= A(r,^)e'\/ rJ -s* 3 * . (5.11) 

Now inverting the Fourier transform in (5. LI), we have 

u(x,i/,f) = j j A{ T ,<;)e-' rtr '™ u \/ rr -^ £ drd<; * (5.12) 

Even though (5.12) gives a perfectly nonreflecting condition at r = L it 
is rather impractical to impose computationally. However, in the pseudo- 
dilFcrenttal operator terminology i\[t 2 — x 2 is nothing but the symbol of the 

pseudo differential operator \J ,y^ — Engquist and Majda proceeded 

with this guideline that the approximations of this symbol i\jr 2 — <; 2 yield 
a family of approximate boundary conditions. In particular, approximations 
around <; = 0 (this has the physical meaning of near normal incidence) gives 
the boundary conditions that we are seeking. For example, the zeroth order 
approximation of i\fr 2 — <; 2 is 

i\/t 2 — <; 2 ir . (5.15) 


Now differentiate (5.12) with respect to x. 

Substitute (5.15) into (5.1 1) and evaluate at x = L to obtain 

« 1) = I J i: 


(5 II) 


(in | 
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The right hand side of (5.15) is nothing but the Fourier transform of -u*. 
Thus (5.15) can be rewritten as 

( 0 d \ 

\ t)x ()t ) 

wliich is the one dimensional condition (1.15). Obviously this will not be 
accurate enough for high resolution to handle two dimensional effects. Fn- 
gquist & Majda proceeded to obtain higher order approximations of the 
symbol i\Jr 2 — < ; 2 . The next order Taylor approximation is 

*\A 2 ~ S 2 ^ ^ . (5.17) 

Substitution of (5.17) in (5.1-1) yields 

u,(L,n,l)=J /■'(’•- l 7 ) /t(f,r) e - r ' +, '» +, ^ n: ^ t rfrdc (5.18) 

Differentiating (5.18) with respect to t we see that 

«i t = Hu ~ ^ liyy , at x = L . (5 19) 

In this way higher order boundary conditions can be generated. One warning 
should be given that the higher order Taylor approximants of the symbol 
i\fr ^ — <;* do not always yield stable boundary conditions. The proof of 
this statement is difficult and will be found in reference 171 . Instead these 
authors proposed Fade’ approximants and found out that they are stable 
for all approximate boundary conditions. The second Taylor approximation 
coincides with the first Dade' approximation and from physical reasoning 
both the boundary conditions (5 10) and (5.11) are stable 

There is another independent work due to Reynolds 20 '. parallel to the 
work of Fngquist and Majda which derives the boundary condition (5 19) 
as a special case. Moreover, this work describes in detail on reducing edge 
reflections. Recently, Keys ' 2 ' 1 proposed a new method again by decomposing 
the wave equation into incoming and outgoing components to obtain a family 
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u =0 , (5.16) 

z-L 



of boundary conditions. In this work, he derived Engquist and Majda’s as 
well as Itoynolds’s conditions as special cases. 

Now suppose cylindrically radiating boundary conditions are imposed. 
Then it is necessary to change the equation (5.5) into cylindrical polar co- 
ordinates: 


n tt — u rr + '7 UQQ -| — a r • (5 20) 

r z r 

This has nonconstant coefficients and needs sufficient modification to apply 
the above theory. Again variable coefficient theory pseudo differential op- 
erators can be used, ft is difficult to summarize this procedure; however, 
there is a similar approach to that which we have given above. First Pade 
approximation of the resulting symbol yields the boundary condition 

/ i) 0 l \ 

\dr dt 2 r) 

to be imposed at sufficiently large distance r. 

ft is interesting to note that this condition can be obtained from separa- 
tion of variables analogous to that of the spherically outgoing solution (5 5) 
In this case we do not have a simple form of D’Alembert’s solution, rather 
it has the form 


(5 21) 


(a o{ e ) + a -M + . 
Eliminating a 0 (fl) in (5.22) we see that 

>t = o(l/r’ ,2 V 


(5.22) 



(5. 25) 


We sec’ that boundary condition (5 21) agrees with the physical one (5.25) 
Eliminating at (0), 0 , 2 ( 0 ) etc., Uayhss Turkel'^ obtained a family of radi- 
ation conditions, which agree only with (5.25) of Engquist Sc Majda Ollier 
higher order conditions differ from each other. Higher order conditions are 
very appealing, but difficult to implement. It will be of interest to both 
applied mathematicians as well as engineers to see how effective are these 



boundary conditions. llariharan and Bayliss^ 10 ! implemented three dimen- 
sional version of (5.23), i.e., (5.1) to a practical problem described below. 
In three dimensions this boundary condition is asymptotically accurate to 



The problem is to solve for sound radiation into atmosphere from a 
cylindrical pipe. There is an incident wave on the left end of the pipe, and 
sound radiates into atmosphere from open end of the pipe (figure 9). 



Z 


Figure 9 Computational plane of sound radiation problem 


The incident pressure wave, in particular, has the form 

p(r,0,p)~f(r,9)e"”*> (5 21) 

where m is called the mode number. The governing wave equations have the 
form 


()n v -I- imw 

— + u 3 + v r = 0 

at r 

+ %.o 

f)t c)z 
<)v Op 

M + 7)r = 0 

Ow im 

_ + — p = 0 

at r 


(5.25) 

(5.26) 

(5.27) 

(5.28) 


These equations are derived from Euler equations of the associated 
flow problem and are available in reference' 101 Due to cylindrical svmmetrv 
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(which one can recognize from these equations) it is sufficient to solve the 
problem in one plane. In this plane as in figure 5.2 the truncated boundaries 
are 71, 7a, and 7 j together with the axis of tne pipe and the inflow boundary 
where a source term is to be imposed. 

There are boundary conditions on the inflow boundary, axis, and on 
the walls of the pipe. But we shall be concerned only with the radiation 
condition on 7i,72,and 75 1'or rest, of the details, the reader is referred to 
reference I lo l again. When the condition ( 5.23) is imposed for p, the acoustic 
pressure, we have 


dp 

fr 





Suppose a point on 7,(1 = 1,2,3) is at a distance R = \/z 2 -f r 2 from the 
origin and tne line from the origin to the point makes an angle a with the 
axis. Then 


dp dp dp . , 

„ = -- cos a + -j- sin a. (5.30) 

dll dz dr 

But from (5.20) and (3.27) we see that 

dp du 

dz dt 

dp du 

f)r = ~dt' 

Thus (5.20) becomes 

-y - ( v cos a 4- v sin a) + = 0 . (5.30) 

This condition, together with other boundarv conditions, was used to 
solve the system (5 25) through (5 28) by a fourth order finite difference 
scheme to obtain solutions reported in figure 5 32. In this figure the vertical 
axis measures the sound pressure levels (dB) and the horizontal avis gives 
angles 0 where the calculations were made at a distance of 10 diameters 
of the pipe. For this situation when the inflow has a time dependence of 
the form VVeiner-IIopf solutions can be computed. We compared our 
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procedure with the work of Savkar and Edelfeltl 17 ! which uses the Wemer- 
( [opf technique. Also, wo compared the soiuttons with some experimental 
data from Vi lie and Silcoxl ,9 L The comparison is shown in figure 5.3 for a 
wave number k = 3 37 and for azimuthal angular dependence of the solution 
e l20 . Weiner- 1 lopf theory and the numerical solutions agree well. However, 
there is a discrepancy with the experimental results, especially near the axis 
for small values of 0. The reason is due to certain uncontrollable factors, 
such as plane wave and lower order mode of e l ° dependence, which arise m 
the experimental situations. The main point of emphasis here is that the 
radiation condition (5.23) is a suitable condition for a practical problem. 



Figure 10 Comparison of results 
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